import math
from sympy import symbols, diff

if __name__ == '__main__':
    # 定义符号和函数
    c, m_b, t_1, t_2, d_s, d_b, h_s, h_b, t = symbols('c m_b t_1 t_2 d_s d_b h_s h_b t')
    f = -m_b*c*((d_b/2.0+2*h_b)/(d_b+2*h_b)*1.0)*((4*h_s)/(math.pi*(d_s**2)*(t_1-t_2)))*t

    # 计算偏导数
    df_m_b = diff(f, m_b)
    df_t_1 = diff(f, t_1)
    df_t_2 = diff(f, t_2)
    df_d_s = diff(f, d_s)
    df_d_b = diff(f, d_b)
    df_h_s = diff(f, h_s)
    df_h_b = diff(f, h_b)

    print("f 对 m 的偏导公式：", df_m_b)
    print("f 对 T_1p 的偏导公式：", df_t_1)
    print("f 对 T_2p 的偏导公式：", df_t_2)
    print("f 对 D_s 的偏导公式：", df_d_s)
    print("f 对 D_b 的偏导公式：", df_d_b)
    print("f 对 h_s 的偏导公式：", df_h_s)
    print("f 对 h_b 的偏导公式：", df_h_b)

    # 获取用户输入
    c = 385
    m_b = 1.137
    t_1 = float(input("输入你测量的 T_1p : (单位: 摄氏度)")) + 273.15
    t_2 = float(input("输入你测量的 T_2p : (单位: 摄氏度)")) + 273.15
    d_s = float(input("输入你测量的 D_s 平均值: (单位: cm)")) * 0.01
    d_b = float(input("输入你测量的 D_b 平均值: (单位: cm)")) * 0.01
    h_s = float(input("输入你测量的 h_s 平均值: (单位: cm)")) * 0.01
    h_b = float(input("输入你测量的 h_b 平均值: (单位: cm)")) * 0.01
    t = float(input("输入你拟合出的冷却速率 T' : (不要遗漏负号)"))

    # 计算各个偏导数的值
    df_m_b_value = -1.27323954473516*c*h_s*t*(0.5*d_b + 2*h_b)/(d_s**2*(d_b + 2*h_b)*(t_1 - t_2))
    df_t_1_value = 1.27323954473516*c*h_s*m_b*t*(0.5*d_b + 2*h_b)/(d_s**2*(d_b + 2*h_b)*(t_1 - t_2)**2)
    df_t_2_value = -1.27323954473516*c*h_s*m_b*t*(0.5*d_b + 2*h_b)/(d_s**2*(d_b + 2*h_b)*(t_1 - t_2)**2)
    df_d_s_value = 2.54647908947033*c*h_s*m_b*t*(0.5*d_b + 2*h_b)/(d_s**3*(d_b + 2*h_b)*(t_1 - t_2))
    df_d_b_value = 1.27323954473516*c*h_s*m_b*t*(0.5*d_b + 2*h_b)/(d_s**2*(d_b + 2*h_b)**2*(t_1 - t_2)) - 0.636619772367581*c*h_s*m_b*t/(d_s**2*(d_b + 2*h_b)*(t_1 - t_2))
    df_h_s_value = -1.27323954473516*c*m_b*t*(0.5*d_b + 2*h_b)/(d_s**2*(d_b + 2*h_b)*(t_1 - t_2))
    df_h_b_value = 2.54647908947033*c*h_s*m_b*t*(0.5*d_b + 2*h_b)/(d_s**2*(d_b + 2*h_b)**2*(t_1 - t_2)) - 2.54647908947033*c*h_s*m_b*t/(d_s**2*(d_b + 2*h_b)*(t_1 - t_2))

    print("f 对 m_b 的偏导数：", df_m_b_value)
    print("f 对 t_1 的偏导数：", df_t_1_value)
    print("f 对 t_2 的偏导数：", df_t_2_value)
    print("f 对 d_s 的偏导数：", df_d_s_value)
    print("f 对 d_b 的偏导数：", df_d_b_value)
    print("f 对 h_s 的偏导数：", df_h_s_value)
    print("f 对 h_b 的偏导数：", df_h_b_value)

    # 计算总不确定度
    u_m_b = 0.00006
    u_t_1 = 0.14434
    u_t_2 = 0.14434
    u_d_s = float(input("输入你测量的 D_s 的不确定度:"))
    u_d_b = float(input("输入你测量的 D_b 的不确定度:"))
    u_h_s = float(input("输入你测量的 h_s 的不确定度:"))
    u_h_b = float(input("输入你测量的 h_b 的不确定度:"))

    u = math.sqrt((df_m_b_value*u_m_b)**2 + (df_t_1_value*u_t_1)**2 + (df_t_2_value*u_t_2)**2 + (df_d_s_value*u_d_s)**2 + (df_d_b_value*u_d_b)**2 + (df_h_s_value*u_h_s)**2 + (df_h_b_value*u_h_b)**2)
    print("总不确定度为：", u)
